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Abstract 

Recent experiments indicate that electromagnetic hysteresis behav- 
ior can be exhibited at the molecular level. Based on these indications, 
a simple dimeric molecule with a hysteresis-like pathway with different 
spatial coordinates for formation and break-up is described and a MD 
simulation using 2-body potentials and switches to form and break bonds 
is implemented to determine whether chemical reaction pathways might 
also exhibit analogous behavior whilst preserving conventional thermo- 
dynamical outcomes. The results of various common thermodynamical 
and kinetic properties are presented, where no unusual thermodynamics 
is observed for the chemical reaction with the loop pathway which might 
be interpreted as not being "'time reversible invariant"' and therefore 
susceptible to manifesting unusual thermodynamical phenomena. This 
system may well model particles that interact via the electro-magnetic 
field to form molecules, since hysteresis behavior is standard and well rep- 
resented in large-scale magnetic and electrical systems, especially in the 
solid state. The potential switching technique circumvents the problem 
of computer intensive three-body calculations which makes this model 
particularly suitable for numerical investigations in both equilibrium and 
nonequilibrium states where the details of a particular reaction mecha- 
nism and its reaction coordinates are not the focus of attention. A new 
algorithm for the conservation of energy and momentum is incorporated in 
regions where the potentials are switched. The thermodynamical param- 
eters determined include the standard free energy, enthalpy and entropy, 
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the activity coefficient ratios, the equilibrium constant and the many en- 
ergy distribution functions of the molecule, all of which are compared to 
the Maxwell distribution function. Both the free dimer and atom particle 
kinetic energy distributions agree fully with Maxwell-Boltzmann statis- 
tics but the distribution for the relative kinetic energy of bonded atoms 
does not, thereby opening to question recent far-from-equilibrium theo- 
ries such as MNET that make use of this presupposition in a fundamen- 
tal sense. The kinetic parameters determined include the rotational and 
translational diffusion coefficients and the Arrhenius parameters. Most 
applications of rotational diffusion seem to presuppose a projection of a 
value about the rotational axis (leading to a cosine dependence) but here 
it is shown that an angular dependence is also a feasible model during a 
first order relaxation process. An NEMD simulation which uses a novel 
difference equation to test for mass conservation is presented which shows 
numerically that the principle of local equilibrium (PLE) is violated and 
is at best only an approximation. The results suggest that although the 
reaction is microscopically loop-like, unlike all the models routinely pro- 
posed, yet the thermodynamics is entirely '"normal"' and yields results 
that do not contradict any of the known laws of thermodynamics. It is 
therefore postulated that reaction dynamics involving hysteresis mecha- 
nisms can occur in nature and may be experimentally verifiable, although 
experimental interpretations tend to construct models that avoid such 
mechanisms. A revision of the concept of "'time reversibility"' to accom- 
modate the above results is suggested. The general design of the reaction 
mechanism also allows for the use of conventional potentials without hys- 
teresis and this will be the object of future investigation. 



1 Introduction 

Recently, experiments have detected the presence of magnetic hysteresis behav- 
ior at the single molecule level ^ E] ; synthesis of such systems are also a hot 
topic of research. Such facts suggest that non-single- valued functions are 
involved in the phase trajectory of the system. A rational extension of this con- 
cept, which has profound theoretical implications is to construct a dynamical 
trajectory where the region of formation of the molecule does not coincide with 
that of its breakdown. There has been a reluctance in the past to consider such 
loop or hysteresis systems because of the absence of experimental evidence of 
hysteresis behavior at the molecular level, and because of the influence of the 
belief of '"time-symmetry"' invariance which discourages such a view, which 
lead to the construction of dynamical pathways which were both single valued 
and which did not have any loop or circular topology; a detailed mathematical 
examination of these common time symmetry presuppositions -so essential to 
physics- has been made^J |S] and it was shown that such views are often not 
warranted or incorrect. This work reports a workable model hysteresis reac- 
tion pathway which leads to thermodynamically consistent behavior, exhibiting 
properties that will require new developments in reaction theory, and it also pre- 
dicts the feasibility of such mechanisms in nature. It suggests a re-definition and 
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extension of the ideas of "'time reversibility"' and "'microscopic reversibility"' 
to cater for the proposed mechanism. Incidentally, the shape of the poten- 
tials and switching mechanism used here is surprisingly similar to experimental 
discussions of the charge neutralization reaction [S] 

k+ + i--»k + i (i) 

except that the discussion does not explicitly mention the crossing over of the KI 
and K + I~ potentials at short distances (high energy) due to the 'time-reversal" ' 
presuppositions referred too above. The existence of a cross-over would make 
the potential mathematically equivalent to the present treatment and there is 
good reason to suppose that such processes can and should occur in electro- 
magnetically induced reaction pathways (such as is manifested in charge-transfer 
and Harpoon mechanisms). The dimeric particle reaction simulated may be 
written 

fci 

2A «=► A 2 (2) 
fc_i 

where k\ is the forward rate constant and k-i is the backward constant. The re- 
action simulation was conducted at a mean temperature which is high, T* = 8.0 
well above the supercritical regime of the LJ fluid. There have been vari- 
ous attempts in modeling chemical reactions with different objectives in mind 
[7| |H1 H3 HOI 1111 IT!?] . Some used generalized models with few details to predict 
the main features experiments might reveal [7] at the reaction coordinate close 
to the transition state (TS), such as what might occur within a solvent-caged 
reaction complex: A-H • ■ • B ^ A • • ■ H-B . This particular pioneering approach 
was further elaborated by Bergsma et al |12| in order to examine the limits of 
validity of TS theory (TST) by not carrying out an ab initio study of all the 
possible reactive trajectories, but by examining trajectories constrained to the 
TS surface because of the limits of computing power. An example of an ab initio 
detailed chemical reaction approach with a 1000 atom system using an assumed 
3 body potential for the exchange process F + F 2 ^ F 2 + F is that of Stillinger 
et al PO] who admits that the procedure is 'very demanding' . At the other 
extreme are generalized studies of hypothetical schemes % such as the 'chemical 
reaction' A + A ^ B + B used to elucidate some kinetic properties. Clearly in 
such models, species A and B must represent complex systems that can be phys- 
ically distinguished; in chemical applications, they might represent for instance 
cis and trans isomers of some compound or they might represent mesoscopic 
species. Some simulations do away altogether with the details of molecular 
dynamics based on dynamical laws [B], replacing them with the Ansatz that 
the details of the interaction between individual particles are not essential in 
the study of the statistical evolution of the system. Such an approach would 
make studies attempting to correlate the details of the dynamics to macroscopic 
properties difficult or obscure, despite the great savings in computer time, and 
therefore does not suite the purposes at hand here. The objectives of the present 
study include: 
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(a) designing a mechanically well defined reaction model with low computational 
demands and where the averaged motions of the dimer may be correlated with 
the macroscopic kinetic and thermodynamical properties and where no anoma- 
lies must be observed in the macroscopic results. Such an outcome would imply 
that the dynamics are reliable enough to be used in other studies 

(b) introducing some degree of complexity to the dimer such as vibrational and 
rotational states for more detailed dynamical investigations 

(c) utilizing the thermodynamically consistent model (as judged by the results 
of an equilibrium simulation) in nonequilibrium simulations 

Here we focus primarily on (a) above. To this end a new general algorithm 
(which will be discussed separately in another planned work) was used to con- 
serve momentum and energy. 

The following essential thermokinetic parameters will be determined and 
discussed in the sections that follow: 

• The thermodynamic equilibrium constant through extrapolating the den- 
sity to zero. 

• The activity coefficient ratio. 

• The standard Gibbs Free energy, Enthalpy and Entropy of the reaction 
through extrapolation. 

• The Arrhenius activation energy and pre-exponential terms, which bears 
no immediate connection to the potential of activation in Fig. ^ an d the 
rate constants of the forward and reverse reactions. 

• The diverse probability distributions for the kinetic energy about the CM 
(center of mass) for all the species, the internal energies of the molecule, 
which do not have a Boltzmann distribution , thereby casting into doubt 
many fundamental theories that assume the opposite (where is is noted 
that the famous Eyring Activated Complex Theory (ACT) does not con- 
sider the vibrational mode for the reaction coordinate to be active, thereby 
not contradicting our results)and these are compared to the Maxwell- 
Boltzmann distribution. 

• Self Diffusion and rotational diffusion constants. 

The method appears very promising for quantitative simulations of real systems, 
and will be utilized in the years ahead for various reaction studies. 

2 The Model 

We examine the dimeric particle reaction given in J2J) above 

2A^ A 2 
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in a range of equilibrium fluid states all well above the LJ supercritical regime. 
This model resembles somewhat that of rcf. except that a harmonic po- 
tential is coupled to the products to form the bond of the dimer whenever the 
internuclear distance reaches the critical value rj between two free atoms A. 
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Figure 1: Potentials used for this work 



In the current study, the potentials as given in Fig. ^ ar e used, but other 
configurations are possible, as verified by direct simulation, such as the excited 
state configuration of Fig. [3 and the reduced distance model model with the 
same spatial coordinates for the forward and reverse reactions in Fig. 01 This is 
a typical reaction potential and it is proposed that a quantitative simulation of 
a simple dissociation reaction of a diatomic gas be attempted. It was found that 
the equilibrium exchange rate of eqn.Elwas very low at lower temperatures and 
changed rapidly at higher temperatures to a saturation level for the latter model 
(Fig. |3J, not making it very suitable for studies where rates of formation and 
breakdown of bonds must be large enough for accurate statistics to be gained 
across the MD cell over a wide range of density and temperature ranges; the 
reason for the slow exchange is in part related to the small reaction or collisional 
cross-section of the molecule. 

The MD mechanism for bond formation and breakup is as follows. The free 
atoms A interact with all other particles (whether A or A2) via a Lennard- 
Jones spline potential and this type of potential has been described in great 
detail elsewhere p^- An atom at a distance r to another particle possesses a 
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Potentials for excited state model 
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Figure 2: Potentials used for the excited molecular state 
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Figure 3: Potentials used for reduced distance molecular model 
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Figure 4: Pressure and temperature distribution across the MD cell 



mutual potential energy ulj where 



ulj = 4e 



for r < r s (3) 



To 



ulj = Oij(r - r c ) 2 + bij(r - r c ) 3 for r s < r < 

ulj = for r > r c 

and where r s = (26/7) i a The molecular cut-off radius r c of the spline 

potential is such that r c — (67/48)r s . The sum of particle diameters is a and 
e is the potential depth for interactions of type A-A (particle-particle) or A-A2 
(particle-molecule) designated (1-1) or (1-2) respectively. The constants and 
bij where given before p|! as 

<Hj = -(24192/3211)e/r2 

bij = - (387072/61009) e/r 3 s (4) 

The potentials for this system is illustrated in Fig.^ Any two unbounded atoms 
interact with the above ulj (1-1) potential up to distance r/ with energy e = 
ULj(rf) when the potential is switched at the cross-over point to the molecular 
potential given by 



u(r) = u mb (r)s(r) + u L j [1 - s(r)} (5) 

for the interaction potential between the bonded particles constituting the molecule 
where u V ib(r) is the vibrational potential given by eq.© below and the switch- 
ing function s(r) has the form given by eq.JJJ) . LJ reduced units are used 
throughout this work unless stated otherwise by setting a and e to unity in the 
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above potential description. The relationship between normal laboratory units, 
that of the MD cell and the LJ units have been extensively tabulated and dis- 
cussed [T3]and will not be repeated here . For the system simulated here with 
the potentials depicted in Fig. the switching function is operative up to 
Tb, the distance at which the molecule ceases to exist, and where the atoms 
which were part of the molecule interact with the (1 — 1) potential Ulj like 
other free atoms;bonded atoms interact with other particles , whether bonded 
or free with the ulj (1-2) potential. The point Tf of formation corresponds to 
the intersection of the harmonic u V ib(r) and ulj curves , and their gradients 
are almost the same at this point; by the Third dynamical law, momentum is 
always conserved during the crossover despite finite changes in the gradient, 
since the sudden change of the force field is between only the two particles 
where the Third Law applies, thereby conserving momentum. Total energy is 
conserved since the curves cross, and errors can only be due to the finite time 
step per cycle in the Verlet leap frog algorithm, which would cause the atoms to 
be defined as molecules at distances r < r/. Similarly at the point of breakup, 
there is a very small (~ 10 -4 LJ units of energy) energy difference between the 
LJ and molecular potentials despite using the switching function in the vicinity 
of the region to smoothen and unify the curves; the small energy differences 
at the cross-over points are less than that due to the normal potential cut-off 
at distance r c where the normal (unsplined) LJ potential is used in MD sim- 
ulations. In order to overcome this problem, a new algorithm (NEWAL) was 
developed, the details of which will be described in another work which conserves 
momentum and energy at these cross-over points. If E p {r) is the inter-particle 
potential (energy) and E m {r) that for the molecule just after the crossover, the 
algorithm promotes the particles to a molecule and rescales the particle veloc- 
ities of only the two atoms forming the bond from Vj to v'i (i = 1,2) where 
v'i = (1 + a)vi + (3 such that energy and momentum is conserved, yielding 
(3 = ~" ( '"^^:^™ 2V2 ' > (for momentum conservation) and energy conservation im- 
plies that a is determined from the quadratic equation a 2 qa + 2qaa — A = 
with a = (vi - v 2 ) 2 ,q= 2 (ml+m 2 ) and A = ( E P ~ Em ^> ' Intercnan S in g 171 and 
p allows for the same equation to be used for break-up of the molecule to free 
particles. For the simulations, success in real solutions for a for each instance 
of molecular formation is 99.9 % and 100% for breakdown-where the A value 
in this instance is very small ( ~ 1.0 x 10~ 4 ). In these simulations, we ignored 
the cases when there was no solution to the quadratic equation, meaning no 
molecules are formed at all, and the interactions are of the (1 — 1) variety. This 
new algorithm coupled with shorter time step (from 0.002* to 0.00005* ensured 
excellent thermostatting, where the thermostating was carried out at the ends 
of the box, as is the case in most real systems. It should be noted that this 
much smaller time scale is not unrealistic as the temperature for this system is 
of the order of 20 — 30 larger than the usual values chosen, and so the trans- 
lational kinetic energy of the particles would scale by the same order. In this 
equilibrium study, the MD cell (which is a rectangular box) is divided into 128 
equal orthogonal layers in the x direction, which is of unit length in cell units. 
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In this method of boundary conditions ^21 j the first 64 layers to the midpoint 
along the x axis are a mirror reflection about the plane parallel to the other two 
axis passing through this x axis mid-point. The y and z directions have length 
1/16 each (cell units). This shape is chosen because non-equilibrium simulations 
will concentrate on imposing thermal and flux gradients along the x— axis, which 
would allow for more accurate sampling of steady state properties about this 
axis |14| . The layers that are mirror reflections about the mid-point plane are 
averaged for steady state thermodynamical properties, leading to effectively 64 
layers. With this algorithm, with only wall thermostatting, we sample each of 
the layers for temperature and pressure changes, and find that the profiles arc 
rather constant, as shown in Fig.^ The heat supply term (per unit time) are 
zero to within the error of fluctuation of energy. Without the algorithm, (but 
with the same time step increment )the center of the effective cell (layer 32 ) 
would have a temperature T* higher than that of the thermostatted end lay- 
ers by over 2 units, and the heat supply term would be significantly negative, 
implying a virtual heating up of the system at the middle due to the potential 
differences due to the switches at the crossover points. The pressure too would 
be unrealistically higher at the center of the cell, which is unphysical. The al- 
gorithm above therefore is very effective in overcoming these problems. Prior 
to this, each layer would be thermostatted to maintain a constant temperature 
and pressure profile, conversion units At regions r < r sw ,s(r) — > 1 imply- 
ing u(r) ~ u v ib(r), i.e. the internal force field is essentially harmonic for the 
molecule and at distances r < r sw ,u{r) ~ ulj, so that the particle approaches 
that of the free LJ type. Concerning the mechanism for the switching, in quan- 
tum mechanical kinetic descriptions, switch mechanisms are frequently used 
for describing potential crossovers 0, but from a classical viewpoint one can 
suggest that the inductive LJ forces due to the particle potential field (with 
particles having a state characterized by state variables slj) causes the internal 
variables at the critical distances and energies mentioned above to switch to 
state Sm when another force field is activated for the atoms of the dimer pair. 
State Sm reverts again to state s^j at distances r^. The shape of the potentials 
and switching mechanism used here is surprisingly similar to discussions of the 
charge neutralization reaction (Hj mentioned in (JIJ K + + I - — > K + I except 
that the discussion does not explicitly mention the crossing over of the KI and 
K + I~ potentials at short distances (high energy), although there is reason to 
suppose that such processes may well occur, since the KI potential curve exists 
at shorter distances well before the crossover point. The following values were 
used here for the potential parameters: 

(a) Current study (Fig. 

Mo = — 10, ro = 1.0, k ~ 2446 (exact value is determined by the other input 
parameters) ,n = 100, Tf = 0.85, r b = 1.20, and r sw = 1.11. 

(b) Excited state model (Fig. [2J| 

Mo = 10, ro = 1.0, A; ~ 2446 (exact value is determined by the other input 
parameters) ,n — 100, tj — 0.85, = 1.30, and r sw — 1.17. 

(c) Reduced distance model (Fig. |3J) 

Mo = —8, ro = 0.6, k ~ 2446 (exact value is determined by the other input 
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parameters) ,n — 100, rj — 0.90, ri> — 0.90, and r sw — 0.90. 

The intramolecular vibrational potential u V ib(r) for a molecule is given by 

u v ib(r) = u + ifc(r - r ) 2 (6) 

A molecule is formed when two colliding free particles have the potential energy 
u{rf) whenever r = rj < ro, at the value indicated in (a) above. This value can 
be defined as the isolated 2-body activation energy of the reaction and has the 
value of 17.5153 at the indicated value of r/. A molecule dissociates to two free 
atoms when the internuclear distance exceeds (which in this case is 1.20). 
The switching function s(r) is defined as 

s(r) = (7) 



where 

s(r) -> 1 if r < r sw 
s(r) ^ for r > r sw 

. The switching function becomes effective when the distance between the atoms 
approach the value r sw (see Fig. Q). 

Some comments concerning the MD potentials are in order. It is generally 
not correct to assume that the potentials in Fig. Q represents the transition 
state theory (TST) potential surfaces; these surfaces can only be derived by 
computing the actual potential of the dimer or free atoms at a known internu- 
clear distance in the presence of all the other species: the zero density limiting 
potentials of Fig. Q cannot cause stable molecules to exist if they were formed 
by excited atoms with total kinetic in excess of the zero density activation en- 
ergy since if energy is conserved, the formed molecule would (except for a finite 
number of kinetic energy values, depending on the model) have to dissociate 
again to the atomic states from which they were formed initially. There must 
be energy interchange at the potential well of the molecular species to remove 
energy so as to prevent dissociation. This is achieved through the presence of 
the temperature reservoir. This reservoir, if it is coupled to the system would 
induce a system behavior whose limit at zero density would not be the same as 
an isolated mechanical system. Likewise, all other state functions of activation 
(free energy, entropy, etc. ) must be computed as functions of all the coordi- 
nates of the particles involved in the interaction. The numerical magnitude of 
these functions cannot be inferred from the isolated potentials above. It is sur- 
mised that these are the potentials that must be used to determine via statistical 
mechanics the various system properties, such as the equilibrium constant and 
the state functions. Here, we extrapolate to zero density at fixed temperature 
to derive these functions, which cannot be inferred from mechanics only. 
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3 Thermodynamic results from equilibrium mix- 
tures 



The reacting mixture considered here were in thermodynamic equilibrium with 
4096 particles. The cell was thermostatted at the ends of the cell maintained 
at the same temperature. Typical runs of 10 million time steps were performed 
per run at each general particle density p (where p is determined as a general 
density irrespective of whether the particle is free or is part of a molecule), 
where the first 200, 000 steps were discarded so that proper equilibration could 
be achieved for our data samples. The sampling methods have been previously 
described ^3] where sampling of all data variables were done each 20" 1 time 
step and where there were 100 dump values where each dump consists typically 
of 5 x 10 5 samples which are averaged. The 100 dump values are then averaged 
again to yield the standard errors of all variables. Dynamical quantities however 
had to be sampled at each time step St* = 0.00005. The thermostatting method 
conserves momentum and registers the energy absorbed at the thermostats |15j . 
All parameters given here are relative to LJ reduced units, sometimes denoted 
by *. 

3.1 Equilibrium constants 

In order to find the thermodynamic equilibrium constant, K eq , the following 
procedure was adopted. The concentration ratio, K c defined as 

K c = (8) 
x A 

was determined as a function of average system density, p where the x's rep- 
resent number density concentrations. For this and other static quantities, the 
temperature was set at T* — 8.0. At very small densities, the system becomes 
an 'ideal' mixture, but as mentioned previously, the limit of the potentials can- 
not be the same as the isolated potentials used in the MD calculations, since if 
this were the case, all the molecules would break up, yielding a net zero value 
for the equilibrium constant at the limit of zero density. As another project, 
it would be of interest to determine the limit at which the equilibrium regime 
breaks down in this thermostatted system, and to elucidate the theory when 
this occurs. There may well be technical difficulties involved in computations of 
very low density systems though . The plot of K c = K c {p) is shown in Fig. 
The accuracy of the K c values varies inversely with p, where in the captions 
sd refers to the number of standard deviations of the standard error. At low 
densities, fluctuations in K c implies that any extrapolative method can be ruled 
out, unlike previously (when NEWAL was not devised) when all the layers were 
individually thermostatted and where a least squares fit n order polynomial ex- 
pansion p(x) = y^_ Q aix 1 to derive the zero density limit of the concentration 
ratio was utilized; the value of n was between 2 to 4. The zero density limit 
Kq where Kq(T*) — K c (p — > 0)is the true equilibrium constant. It is clear that 
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in this system K and K c in general differ significantly; it serves as a warning 
that in general, one cannot ignore activity coefficients in the calculation of such 
properties in model systems and theoretical demonstrations. In the present 
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Figure 5: Variation of concentration ratio K c with p, the system number density 
at LJ temperature T* = 8.0 with sd — 3 at p — .03 and sd = 50 at p = .9 

study, it was discovered that at very low densities, fluctuations are significant 
as shown in Fig. © for the case of a run at T* — 8.0. The method used in the 
present case is to take the mean value of K c for very low p values ranging from 
0.03 to 0.09, for about 12 values at any one temperature and to approximate this 
as Ko(T*). The fluctuations show that in this range of density, the system has 
" 'saturated" ' itself in that all the p values yield approximately the same mean 
K c . The results derived for T* = 8.0 are 

K eq = lim K c = 0.0610 ± .002 LJ units. (9) 

In previous studies prior to NEWAL implementation, using polynomial extrap- 
olation, a value of 0.050 ± .001 was derived. Knowing this value, we calculate 
the activity coefficient ratio, <E>, for the other densities at the same temperature 
by using 

K eq = K c ^ - K c $ (10) 

7a 

The ratio of activity coefficients $ is shown as a function of density in Fig. Q . 

It is clear from the $ ratio that for normal densities, the equilibrium reaction 
mixture is highly non-ideal, which may be expected due to the large differences 
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in the LJ energy well for the molecule and the atom (see Fig.Q. It is probably 
a poor approximation to use ideal models for test systems in reactor design, 
which is often the case. Further, the above technique allows for the general 
determination of activity coefficient ratios via simulation. The determination 
of separate activity coefficients is a challenge. One real problem is the fact 
that molecules, in the equilibrium state cannot exist in isolation. In mixtures, 
either the reaction goes to completion, or they do not react in simple theory of 
mixtures. In these cases one might postulate separate ideal states for the "pure" ' 
components, but in the present elementary case, for any one temperature, there 
is a finite value for Kq meaning the presence of all components in a system 
at equilibrium. It it therefore a challenge to find a suitable model or concept 
to solve this problem with cycle changes. Even if a hypothetical state were 
defined, one must still design the route or cycle taken to the equilibrium state 
which consists of product and reactant species. The derivation would require 
a scries of very elaborate and detailed computations and is not attempted here 
since it is not immediately relevant. The rate constant is a defined quantity, 
which accords with the standard form below. The overall rate of reaction r may 
be written in terms of the experimentally determined forward rate (r\ = k\x\) 

for the process 2 A — \ A 2 and backward rate (r_i = kiXA 2 ) for the process 

A2 — > 2A as r — r± — r_i = k\x\ — k-\XA 2 At equilibrium r = 0, and so 



xa 2 _ _ki_ 

The ratio of rate coefficients is the concent rat ion ratio K r where 



in) 



*• - h <12) 

To verify the above equilibrium constant independently from kinetic measure- 
ments, we can extrapolate to zero density p the values for r\jx\ = Q = k\ and 
r_i/iA 2 = R = fc_iThe rates were calculated independently from the program 
by monitoring the number of bonds formed or broken for each time step St* and 
averaging this quantity over the 10M time steps. Then the relevant equations 
are 

'Q\ T ^ limQ(p^O) Q° 



The plots of Q and R at low densities are given in Fig.lJSJ. 

As for the direct determination of the equilibrium constant, fluctuations 
imply an averaging at very low densities to derive the limits. The results with 
the estimated errors are 

lim Q = Q° = 0.870 ± .006 L.J. units (14) 

lim R = R° = 14.32 ± .1 L.J. units. (15) 
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Figure 8: Low density values of Q and R variables. 



It will be noticed that at very low densities, we would expect the errors due to 
the breakdown process to be very much higher than that due to the formation 
process since the number of dimers tends to a very low number. The ratio of 
these values gives the the true equilibrium constant directly from kinetics as 

fci 

K eq (kinetic) = lim — - = 0.061 ± .001 L.J. units (16) 

An excellent agreement with the results from the previous equilibrium analysis is 
found, where the method used for the determination of the equilibrium constant 
differs. This agreement indicates that the system is in a steady (equilibrium) 
state and that the simulation method is fairly coherent. The Q and R functions 
at other densities are given in Fig.©. 



3.2 Standard states 

We use the form AG°(T) = — kT In K eq to determine the standard free energy 
state AG°(T) of the dimer reaction. The justification is that we can choose the 
standard state to be at constant pressure (of zero value) for the standard state, 
so that the chemical potential standard state for each species is only a function 
of temperature, so that AG°(T) is strictly only a function of temperature |161 
p. 177-179]. We repeat the same process as described above in section (I3.1|) for 
T* = 8 for different temperatures (from T* = 4 — 20. Each determination 
required at least 8 runs at varying low densities. It was found that at low 
temperatures, the fluctuations were greater, as shown in Fig (1 1011 where the 
variation of K eq versus 1/T is given. The curve used to determine the other 
standard state functions was the Gibbs free energy curve , given in Fig. l(TT)l. 
For this curve, the error bars (except for the first data set) all refer to the 
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Figure 9: Variation of Q and R variables with density. 
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Figure 10: Variation of equilibrium constant K eq with 1/T for fixed average 
system p = 0.70 with best fit line 
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errors relative to the least squares fit of a quadratic curve to the simulation 
result. The fit is rather good. The standard entropy A5°(T)is derived from the 
thermodynamical entity [161 eqn. 6.34, p. 182] 

Clearly to use (|17l) . we must know AG°(T)as a function of temperature T. We 
write therefore a simple quadratic equation with p coefficients as follows 

AG°(T) = p(l)T 2 + p(2)T + p(3) (18) 

The non-linear least squares method yields p(\) = —0.0233441, p(2) — 
1.0531305, p(3) = 15.46544989 with an overall uncertainly of the free energy as 
approximately ±0.3. Differentiating 1)18(1 yields the entropy as AS = — (2p(l)T+ 
p(2)), which is linear. The standard enthalpy AH is given at constant temper- 
ature by the entity ^3 p. 183] 

AH°(T) = AG - TAS° (19) 

which therefore means that the standard enthalpy is given by AH° = —p(l)T 2 + 
p(3). It can be verified that this expression and that for AS recovers the 
quadratic (fH^I . 
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Figure 11: Variation of the standard Gibbs Free Energy AG with temperature. 



The plots for the standard entropy and enthalpy as functions of temperature 
are given in Fig. (fT2*|l . 
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Figure 12: Plot of standard enthalpy AH°(T) and entropy AS°(T) from 
AG°(T) quadratic curve fit with the temperature T*. 



The essential point here is that the standard entropy is negative, as it must be 
at moderate to low temperatures since the free particle state has a larger phase 
space than the corresponding dimer. It may appear counter-intuitive that the 
standard enthalpy is positive. It must be pointed out that at these temperatures, 
the particles are not trapped at the bottom of the potential well, and that the 
activation energy is positive, and that the internal potential energy at the point 
of formation of the molecule is not lost, but is converted to internal kinetic 
energy, leading to the break-up of the molecule. A quantitative treatment of 
these terms has been attempted ^7]. It must be concluded that the simulations 
are able to determine the standard states without having to construct extremely 
detailed cycle diagrams; further, the simulation can also check on the correctness 
of the cycle diagrams used to determine standard state values. 

3.3 Activation energies 

From the way the algorithm was constructed for molecular formation, the molec- 
ularity of the elementary reaction is 2 leading to a single second-order reaction 
of formation, and for the dissociation of A2,a first-order reaction results since 
the molecule can only exchange kinetic energy with all other particles within the 
system without further reactions to the dissociation limit. A frequently used 
model for the kinetic constant k% for these rates is due to Arrhenius, which has 
the form 



where the rate constant is a function of the temperature only and where Ai is 
ideally not temperature dependent. It should be noted that the Arrhenius equa- 




(20) 
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tion is strictly valid for 2-dimensional systems where the pre-exponential factor 
is independent of temperature and where the exponential factor exp (— ^) rep- 
resents the fraction of molecules having energy in excess of E$ JSji where Ei 
is usually understood to be the activation energy. The reason why this form is 
so durable is that the exponential term represents the fraction of excited state 
atoms, and this term dominates over the pre-exponential term with temperature 
variation, which give the impression of constant Ai factor for the plots. The 
rate constants for the forward fciand reverse reaction were plotted versus 
1/T for the given density of p = 0.7 and was found to be reasonably linear 
(Figs. (|13I14|1 . with the activation energies for the forward and the backward 
reaction rates (Ei and E-\ respectively) and the corresponding collision factors 
(Ai,A_i) determined approximately as 
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Figure 13: Variation of natural logarithm of forward (product forming) rate 
constant k± with reciprocal of temperature for p = 0.7 



E x = 21.40 ± .10 LJ units, A 1 = 3.50 ± .2 LJ units 

E-i = 7.26 ± .02 LJ units, A_ a = 2.70 ± .04 LJ units 

There are two separate rate constants here, for first and second order. The 
second order forward rate constant k\ has a form given by 

/8kT\ 1 ^ 2 ( e* \ / e* \ 

J =A iaq >{--) (21) 
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Figure 14: Variation of natural logarithm of backward (product disintegration) 
rate constant k-i with reciprocal of temperature for p = 0.7 



Very roughly, if the mean temperature for the plot (which spans from 4 to 20) 
is 12, then l|21|) above yields for the given value of A\ b max = 0.9153.. which is 
reasonably close to 0.85, the theoretical value. However, e* = 21.40, which is 
higher than 17.5153, which is the set simulation potential value for the formation 
of a molecule. Since we can expect a yet greater accuracy for the determination 
of e* as compared to Ai due to the domination of the exponential terms, it may 
be safe to suppose that other factors contribute to the true activation energy 
other than what is described by simple collision theory (SCT). Future work will 
attempt to determine what other energy factors are implicated in e*; currently, 
SCT views this energy as a pure mechanical work energy, which obtains at the 
molecular level. Similarly, variation of Ai with various energy terms cannot be 
immediately ruled out. Generally, the above values do not bear a direct re- 
lationship to the isolated 2-body potentials of Fig. JIJ, but nevertheless some 
approximate correlations are evident; E% is somewhat close to the isolated ac- 
tivation energy 17.5153 measured from the free atomic states, and likewise E-\ 
is somewhat close to the energy difference from the bottom of the molecular 
potential at —10 to the potential at n, which is approximately —1, a distance 
of approximately —9 energy units. However, for a first-order reaction, a differ- 
ent interpretation for energy differences obtain than from that due to SCT for 
instance, which is concerned with bimolecular processes; the first order interpre- 
tation is that the molecule decomposes when it overcomes an energy activation 
threshold, and the fraction of such molecules is reflected in the exponential term, 
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the pre-exponcntial term reflecting the mechanism of the decomposition. 

4 Results from equilibrium dynamical trajectory 
analysis 

This section concentrates on variables which had to be sampled at each time 
step of duration St = 0.00005* in order to compute the property of interest:the 
rate of reaction in the previous section above is also based on instantaneous 
sampling but more properly belongs to topics associated with equilibrium. Of 
importance in noncquilibrium and kinetic studies are the values of the diffu- 
sional coefficients, reaction correlation coefficients and the energy probability 
distributions, where if the principle of local equilibrium (PLE) obtains imply 
that we may approximate the values computed in an equilibrium simulation 
for those in a nonequilibrium volume element having the same state variables. 
Examples of these quantities (which can also gauge the appropriateness of the 
model for nonequilibrium studies) are provided. 

4.1 Rotational diffusion constants 

Although connected in some ways to diffusion, a somewhat unconventional 're- 
orientation' diffusion function (cos 4>(t)) has been defined [7] where <j)(t) is the 
angle between R(0), the unit internuclear distance vector of the dimer at t = 0, 
and R(t), the same unit vector at time t. Such a definition might have appli- 
cations in conjunction with their being part of transform functions eqs.(17)- 
(20), p. 211], where the postulated exponential decay of this function when acting 
as a kernel of the transform could force convergence of the function being con- 
voluted. It is found that the exponential decay assumption in cos(<^(i))is a fair 
but not perfect fit, perhaps implying that another type of theory for "'rota- 
tional diffusion" ' constants may yield even better fits with the experimental 
curves. We provide one such example (arccos(t)), an approximation to (0(t)), 
which provides a far better fit and therefore is a candidate for a stochastic the- 
ory of rotational diffusion. This then is another area for research. It must be 
mentioned, however, that the theory of '"rotational diffusion"' as developed by 
P. Debye and others ^3 p.81-84,esp eqs. 49] etc. makes use of '"dissipation 
kinetics"' where a constant torque M is balanced by a inner frictional force C 
parameter, so that M = £4| , where 9 is an angular displacement. Such a 
theory leads to a relaxation in the distribution function / by a factor ip(t) given 
by tp(t) = exp— =^-t so that for a particular orientation angle 6, f has the 
form / = A [1 + Cip(t) cos 6]. The mean dipole moment of the entire sample 
also decays with the same rate as with t/j. It is not immediately clear that the 
orientation angle must also relax according to a first order rate law. If the ef- 
fect is a projection of an orientation onto an axis,then this would correspond to 
the result given by Allen et al (op cit). O'Konski and Haltner [20] have char- 
acterized TMV (virus) by studying the birefringence relaxation rate written 
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5 = S cxp(— t/r) where t is the initial value of birefringence [201 eqn 3, p. 3607]] 
and the '"rotational diffusion coefficient" ' Dh is defined here as Dh = l/6r with 
an additional factor of f/3 to that of Allen. Most of these theories supposes 
that even at the molecular level, one can use frictional coefficients as for macro- 
scopic systems where the retarding force is linearly proportional to some form of 
velocity of the system , the constant of proportionality involving the frictional 
coefficients [21] ■ More recent studies experimental studies of rotational diffusion 
|22l I23*| assume a first order relaxation of fluorescent directed intensities of the 
chromophore of the molecule with the rotational diffusion constant defined as 
in "20I- To show that the results obtained is typical, we graph the functions as 
defined by Allen et. al [7]. The method used here to determine {cos <j)(t)) is to 
create a table whenever a molecule is formed which maps out for each increment 
in the time step i the value of cos </>(«) until it disintegrates: for each i th time 
step there exists for each sampling subinterval M (M being a variable) values of 
4>{i) due to other molecules which have existed, and the average value for each 
sub-interval is computed as (cos^(i)) = J2jLi cos <j>j(i)/M, According to Allen 
et al (op cit), the function decays as 

(cos </>(*)) = Aexp(-t/n) (A = 1) 
with linearized form 

In ((cos 0(f)}) =-t/n (22) 

where the '"rotational diffusion"' coefficient D r is given by D r = The 
results of the simulation is graphed in Figs. I|f 511 7fl . Fig. "TBI 

graphs the proposal found in 0. ft is clear that there is an initial chaotic 
regime, followed by a very slow decay of approximate form Aexp(t/r r ),(A — 1 
if we measure the time from the end of the chaotic regime onwards; fitting this 
portion of the curve from the 400 — 800th time step to the above exponential 
yields r r = 1.38 ± .02LJ units. A 'rotational diffusion constant' D r — may 
be defined and the value obtained is D r = 0.36 ± 0.01LJ units. The shape of 
the (cos 0(i))curve resembles that described in [7j (where the 'initial chaotic 
region' is mentioned) implying a somewhat typical rotational motion, but it is 
clear from the figure that even in the fitting region, there is an apparent concave 
shape, as the tangent line makes clear. Nevertheless, for the sake of parametriza- 
tion, this particular definition is used to derive the diffusion constant D r data 
at other regimes of varying p (at constant temperature) in Fig. (|17(l and for 
varying temperature (at constant p) as depicted in Fig. I|16|l . In these figures, 
the same method of determining D r was used as for the above determination 
of D r at p = 0.7 and T* — 8. As with the case of rectilinear diffusion mo- 
tion D t = BkT, where B is the density dependent mobility coefficient, which 
is the steady state velocity acquired per unit external force [""""fl sec.l4.4,eqns 
(2)-(ll). p. 464-465], we obtain at fixed density p a linear relationship with tem- 
perature, suggesting a similarity or isomorphous theoretical construct in relation 
to rotational motion. Noting that different thermodynamical variable regimes 
are associated with different error margins when determined experimentally, we 
also notice an approximate linear correlation with density at fixed temperature. 
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Figure 15: Variation of natural logarithm of cos <j) orientation function with time 
at T* = 8.0 and p = 0.7 



From the rectilinear equation, this would be the case if the mobility coefficient 
B were inversely linearly related to the density of the medium, which is a very 
reasonable assumption at higher densities (p* — 0.75 — 1.0). The figures show 
that the change of the diffusion constant with p at fixed temperature is much 
less dramatic than with temperature at fixed p. 

Fig. (|18f) gives a clear indication that the long-time correlation is linear 
concerning time and the logarithm of 8, and so one can also derive a rotational 
diffusion theory where not a projected value, but rather the actual angular 
distance relaxation is a first order process. This, at any rate is what the model 
here depicts. 

4.2 Self-diffusion coefficients 

In these simulations, the mean lifetime of the molecules vary in the region of 
24,000 to 2400 time steps as the corresponding temperature varies from T = 
4.0 to T = 8.0. The accurate determination of the three dimensional (3-D) 
self diffusion coefficient D s for any particle requires the determination of the 
integral of the long time limit of the velocity autocorrelation function, or the 
equivalent Einstein expression of the mean square displacement at infinite time 
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Figure 16: Variation of D r with temperature at constant p = 0.7 
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Figure 17: Variation of D r with density at constant temperature T* 
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Figure 18: Variation of natural logarithm of arccos(cos <j>) orientation function 
with time at T* = 8.0 and p = 0.7 



with respective forms 

1 [°° 

D s = -J dfc(v,(t)-Vi(0)) (23) 
and 

2tD s = l^| ri (i)-r,(0)| 2 ) (i-foo) (24) 

respectively. We overcome the infinite time problem here by determining the 
diffusion coefficient according to (|24|l at the time of breakup tb r ,i of molecule 
i (where the time is when the molecule is formed), thus allowing for the 
maximum time possible before -D s ,m,i is computed (where m refers to the 
dimer). Likewise, we can monitor the time spent as a free particle of any labeled 
atomic species (j), and determine the self diffusion coefficient -D s ,a.j (where a 
refers to the atomic state). The molecular self diffusion coefficient is the aver- 
age of all molecules determined during the dump interval, and lastly the 100 
dump values for the entire run is averaged to provide an estimate of uncertainty. 
Similarly, a labeled particle is used to determine the atomic diffusion coefficient 
based on the time spent as a free, non-bonded particle. The results for this 
supercritical fluid are given in Figs. (|19I20I) . The curves in Fig. (1 19(1 appear 
very linear, verifying the formula D s = BkT, according to previously developed 
theories ([21], eq.(49) ) especially at lower temperatures. The ratio of molecular 
to atomic diffusion constant is relatively close to 0.50 everywhere. The mass 
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Figure 19: Self diffusion coefficients at varying T* and fixed p = 0.7, where A-A 
refers to the dimer and A to the atom, and D denotes the self diffusion coefficient. 
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Figure 20: Diffusion coefficients at varying pand fixed T* = 8.0 
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of the molecule is twice that of the atom and approximately twice the diame- 
ter, leading to this approximate ratio. The actual theoretical prediction due to 
size, energy interaction and mass effects is not well developed, and no exten- 
sive data are available for even non-reacting systems. The reactive system here 
depicts values of the diffusion coefficient which is does not differ significantly 
for systems which do not react. In one study |2H1 P-2044 Table V]of solute 
diffusion in a solvent, where interactions are solvent-solvent (1-1) and solvent- 
solute (1-2) only, (i.e. no (2-2) interactions) the L2 system has the following 
Lennard- Jones parameters — = 2; ^ — 4: ^ = 2 leading to the diffusion 

1 mi 7 en ' <7ii 

coefficients D\ = 0.063 and L>2 = 0.017 (accuracy not specified) and for the 
S2 system, the Lennard- Jones parameters ^ — — \'- §77 = \ lead to 

the diffusion coefficients D\ = 0.082 and D2 = 0.190. for the same mass ra- 
tio, the diffusion constant ratios vary from 0.27 to 0.43 for very different and 
extreme e — a combinations where the variation with temperature is not signifi- 
cant for these ratios based on the scanty information of the graphs drawn; here 
e = 1, andfT = 1 throughout. These ratios are not too different from the ones 
reported here. The variation of the diffusion constant with density is much less 
dramatic than for the temperature according to Fig. I|20(l with a slight decline in 
diffusion constants with increasing density, as is to be expected as the mobility 
would decrease. The errors appear large because the variation of the coefficients 
with varying density is relatively slight for fixed temperature. 

4.3 Energy distribution histograms and non equilibrium 
results 

It is of interest to compare the theoretical Maxwell distribution of the species 
to the distributions derived from simulation since fundamental deductions can 
be made. We also produce more results for a non-equilibrium simulation with a 
novel difference equation which can be used to check for conservation of matter 
to determine whether the principle of local equilibrium is indeed a principle 
or merely a very good approximation for describing general thermodynamical 
systems (whether reversible or not). 

4.3.1 Probability histograms 

These are provided in Figs. (j2H25|) for the translational kinetic energies of the 
different species probed as well as the total internal energy of the dimer. These 
distributions are plotted together with the Maxwell distribution relative to the 
apparent temperature determined from (|27|l . The comparisons provide clues to 
the following: 

• Shape of the probability function P could perhaps be used to determine 
whether the assumptions used in theories is reasonable or not. The shape 
even for this equilibrium system is not always Gaussian, and so there is 
no reason to assume a priori that non equilibrium systems must conform 
to a Gaussian distribution where certain internal variable are concerned. 
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• Providing a rationale for extending the theory of equipartition in an equi- 
librium system where the temperature relative to a particular kinetic en- 
ergy coordinate is not the same as for the total system temperature de- 
termined from standard equipartition. Such a possibility seems to be 
supported by the evidence below. 

For a given Hamiltonian Ti weakly coupled to a heat bath where 

JV 2 

W = X;^- + V(n,r3,...r B ) (25) 

i—l 

where V is the position variable r dependent potential, the probability den- 
sity function per unit area of phase space (p, q) is 

exp-/3H 

P(p,q) = — 7? — (26) 

where the partition function Z has the form 

/ e -/m dpdr 

2 = m ' 

The separability of the Hamiltonian above for the momentum p and position 
variables r which is of the same form as our chemical system Hamiltonian leads 
for large N to the exact result (in 3 dimensional systems) (usual laboratory 
units) 

N{^-)=yyj{2m i ) (27) 



which is the method used to determine the system temperature here. The 
momentum coordinates Pi refer to all atomic species, whether bonded or not. 
The Gibbs postulate can be directly tested for the chemical reaction system. If 
this postulate is valid for loop-like hysteresis systems, then the time trajectory of 
any indexed particle I must also yield, when averaged over a very long time the 
result (in 3-D) SkT/2 = pj/(2mi) whether the particle is bonded or not over the 
trajectory equally weighted for all the states that it traverses. Integrating the V 
function in l|26(l above over all equal energy values, the Maxwellian probability 
density function results, and is given per unit energy increment by 

P = 2^ifcTy /2 e 1 /2 exp _ ( -±. ) (28) 

Ea, 128|l is the standard form used for the absolute velocity distribution func- 
tion since the energy e ot d 2 for velocity v. An apparent temperature param- 
eter < T >x is computed here for some species X and is defined such that 
3<1 2 >x = ^2mr) wnere m x is the mass of species X and px is its momen- 
tum variable. This parameter is clearly not well defined as a temperature if 
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it does not obey the equipartition result above for the obvious reasons con- 
nected to conjugate transforms. In statistical thermodynamics, the total system 
Hamiltonian H = Y^i=\Vl '/^ m + 2~^i<j V{ r i — r j) leads to the density-in-phase 
having form p(p, q) cx exp[— H(p, q)/kT and so for systems with separable 
coordinates, each kinetic energy coordinate E^i = pf/2m and potential form 
V(\ri — rj\) will have the above Boltzmann distribution. However, the "in- 
ternal coordinates" during a chemical reaction or other process refer for ex- 
ample to an artificial aggregation such as the center of mass (CM.) velocity 
and position for particles k, I forming a molecule which is not permanent e.g. 
Pj = Pfc +Pi (k 7^ Rj = m + m ( r *i + r z) need not have Boltzmannized dis- 
tributions. Permanent aggregated states can be expressed in terms of canonical 
transformations Q = Q(p,q), P = (p, q) 27 Chap. VII] and the new Hamilto- 
nian that results must by ensemble theory be subjected to the density distri- 
bution described above. But for systems which are described by '"internal"' 
coordinates of a non-permanent nature ( in the sense that the forces between 
the particles cease when the molecule decomposes) and which does not refer to 
the system Hamiltonian, no general theory exists, and no presuppositions can 
be made to regarding its density distribution. Nevertheless, theories purport- 
ing to be fundamental have been created that assumes the Gaussian density 
for internal variables to be true [281 12"3] without clear qualification concerning 
the situation when this condition obtains. A clear-cut counterexample will be 
provided which therefore opens to question the aforementioned theory. Further- 
more, the principle of local equilibrium has been proposed as essential for 
these new theories, and another counter-example to this is also provided, this 
time from a non-equilibrium simulation. In other words, basic simulation is able 
to determine the veracity of theories, and in particular, the hysteresis system 
described here does not support the novel theoretical developments in " 'meso- 
scopic"' level thermodynamics. The total internal energy coordinate (TIEC) 
and the internal kinetic energy coordinate (IKE) are not Gaussian distributions 
for equilibrium systems according to the simulation result discussed below. Of 
great theoretical interest is that for cases of non-permanent coordinates, some 
types of distributions are essentially Bolztmannized, others are not. It would be 
of great significance and interest to provide criteria which can predict when a 
Boltzmann distribution can be expected. The apparent temperature parameter 
< T >x may well qualify as a temperature in an extended equipartition scheme 
if there is agreement with the Maxwellian distribution even if this temperature 
does not correspond to the unique system temperature < T > sys . Here the de- 
gree of agreement with the Maxwell distribution is either very good (in some 
cases) , or rather bad. It would be of great theoretical interest if some form of 
relationship between the apparent temperatures could be made on the basis of 
internal energetics. The uncertainly (unless stated otherwise) is of the order as 
given in the error bars of Fig. (|25|l which is at 100 standard error units and which 
would not feature in any figure where errors are typically quoted at 3 standard 
error units. This figure corresponds to the TIEC distribution. The errors in 
the temperature are are given in Figs. (|21l27|l . Fig. l|21|) shows that the center 
of mass (CM.) kinetic energy follows quite accurately a Maxwellian P function 
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Figure 21: P functions for translational kinetic energy of A2 about center of 
mass at system temperature set at T* — 8.0 and p = 0.7, with apparent temper- 
ature of molecule indicated. 



with a temperature parameter higher (T* = 8.33 rather than T*8.0 ) than the 
system temperature. The fact that the shape is Maxwellian at the indicated 
temperature parameter does seem to imply that theories may be be developed 
within an equilibrium system with different coexisting temperatures provided 
that these parameters require that a Maxwellian form regarding shape prevails, 
and after that stage one perhaps might also be able to propose generalizations to 
temperature not requiring a Maxwellian distribution; but a proper theory would 
have to begin from first principles which can subsume without contradiction the 
previous axiomatics, including the Zeroth Law. Another inference is that the 
temperatures have definite values (or limits), since the degree of scattering is 
relatively low;hencc one might expect some type of stochastic averaging which 
yields exact values (limits). The other important scientific question is the expla- 
nation of the shift of " 'temperature" ' < T >x for such Boltzmann distributions 
for non-permanent aggregates. An atom bonded to a molecule does not have a 
clear Maxwellian shape, as is evident from Figs. i|22l23[) since there is interference 
from the internuclear potentials. The graph in Fig. (j22|) computes the absolute 
kinetic energy (K.E.(l))of the particle with respect to the MD cell or AKE , 
whereas Fig. (|23(l refers to half the relative kinetic energy and half the trans- 
lational kinetic energy about the CM. of the bonded pair, where the relative 
kinetic energy ek.e.rei., is written as ek.e.rei. = ^(^1 — i"2) 2 = h^ 2 for any two 
bonded atoms 1 and 2, where the reduced mass a is given as - = — + — and 

■ " fj, mi 7712 
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Figure 22: P functions for kinetic energy of any atom A bonded to A2 at system 
temperature set at T* = 8.0 and p = 0.7 with system temperature indicated. 



where the intermolecular axis vector is r = ri — r2- The total internal kinetic 
energy IKE is also defined as the relative kinetic energy of a bonded pair, given 
as ek.e.rei. as above. The AKE averages ||(vi — v 2 ) 2 = \ {vi 2 + v 2 2 — vi.v 2 } 
whereas the kinetic energy about the C.M. ((KCM) averages the expression 

|.2 ^ V1 "1~7 = i {vi 2 + v 2 2 + V1.V2}. Adding these expressions and then di- 
viding by 2 would lead to convergence of the result to that for AKE, which is 
what is presented in Fig. (|23|1 as K.E.(2), which is almost the same graph as for 
Fig. (|22Jl .The reason for this computation was to check for consistency of result 
for the two different sampling techniques. 

The IKE distribution, that of an internal coordinate, is clearly non-Gaussian, 
as depicted in Fig. (|2"l}> . This result is not consistent with the assumptions of 
mesoscopic non equilibrium thermodynamics. |281 129| . 

TIEC defined above refers essentially to the vibrational and rotational kinetic 
energy of the molecule E tiec since the translational kinetic energy about the 
C.M. has been factored away where 

E tiee = V(\r l -T J \) + !^- (29) 

where V^dri — r j | ) = uq + -^k(r — ro) 2 . Hence the intermolecular potential would 
play an important part in determining the motion along the intcrnuclear axis, 
with the environmental potential due to other particles playing a moderating 
role by introducing stochasticity to an otherwise plainly mechanical system. 
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Figure 23: P functions for average kinetic energy of atom A by K.E.(2) method 
at system temperature set at T* — 8.0 and p = 0.7 with total system tempera- 
ture indicated. 




Figure 24: P functions for total internal kinetic energy (IKE) of the two 
bonded A atoms about the internuclear axis at system temperature set at 
T* = 8.0 and p = 0.7 with system temperature as indicated. 
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Figure 25: The total internal energy coordinate (TIEC) distribution as given in 
the text. The error bars are for 100 standard error units. 



The probability of occurrence of a state is proportional to the time spent at any 
configuration, and with a harmonic potential, most of the time spent will be at 
the turning points in simple harmonic motionfin the molecular potential used 
there is a 'dissociation hump' just prior to the dissociation limit, leading to a 
departure from the Maxwell distribution;other reasons for departure form the 
distribution include the dissociation itself, precluding higher energy states from 
being accessed. It is clear that the distribution in Fig. I|25|l is non-Maxwellian 
and accords well with the shape of molecular potential energy function , with its 
humped potential near the distance of dissociation. This model has been used as 
a classic description of equipartition. If the particles were bonded permanently, 
this quantity would have a canonical distribution, which it clearly does not 
because bonds are formed and broken at a rate that precludes adjustment to a 
Gaussian probability factor. This distribution , which also refers to an internal 
coordinate for total internal molecular energy, is not consistent with some recent 
non-equilibrium theories |28l I29| . 

Noting that the accuracy of the single particle is reduced by a factor of 
w 4000 (the number of particles in this simulation), we find that the Gibbs 
postulate seems to be verified in terms of the shape of the P function (which 
appears Maxwellian) as well as the computed value of the temperature with the 
error estimated as ±0.1 by studying an atom of fixed label (no. 29) as it forms 
and breaks bonds with neighboring molecules, as shown in Fig. (|26l) Clearly the 
time average of dynamical properties for this particle would equal the ensemble 
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Figure 26: P functions for kinetic energy of fixed indexed atom A which 
cither is bonded to some A2 dimer or not at system temperature set at 
T* = 8.0 and p = 0.7 with apparent temperature of atom indicated. The un- 
certainty here is 3 standard error units. 
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Figure 27: P functions for kinetic energy of free (unbonded) random atom at 
system temperature set at T* — 8.0 and p = 0.7 with apparent temperature of 
random atom indicated 



average. We notice that the reduced accuracy of the sampling is reflected in the 
greater scatter of the P function points. 

Finallysince the molecular P function has been mentioned, it would be inter- 
esting to compare it to the case of a random, but always free A particle which 
is given in Fig. (|27fl . where the determined temperature is slightly lower, (to 
within the error limits) than the system temperature, and where the shape of 
the P curve is Maxwellian. This particular species type cannot fulfill the Gibbs 
postulate because its trajectory is confined to those areas where there is no 
molecular formation, and so its time averaged properties like the temperature 
need not necessarily equal that for the system as a whole as determined from 
the equipartition principle. We can conclude that the energy subsystems that 
can be chosen for devising a theory of unequal temperature distributions in an 
equilibrium system which all have a Maxwellian probability profile include at 
least the following candidates: 

• Translational k.e. about CM. for A2 

• Fixed indexed k.e. of particle A (in both free and bonded state) 

• Random, always unbonded k.e. of particle A 

The following is suggested as a result the above observations. 
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Conjecture 1 // the random forces are external to the system, and they all 
have the same force law when acting on the particles of the system which may be 
different from the force law for internal forces acting on the particles of the same 
system, then the kinetic energy of the CM. would have a probability distribution 
that is Maxwellian. 

The above conjecture is weak and must be strengthened by a more rational 
theoretical approach using stochastic calculus. 

4.3.2 NEMD results 

Figs. (|28I29|) are the flux and divergence of the flux for "'Case 2"' simulation 
where a temperature gradient across the MD cell is imposed together with the 
making and breaking of bonds at the ends of the cell leading to a molecular flux 
according to the thermodynamical conditions and details given in |14| . The cell 
is broken up into 64 layers along the X-direction and the thermostats are placed 
at the ends of the layers. Fig. I|28|l has overlapping error bars with magnitudes 
that do not change significantly over the range where the fluxes are evident. 
The stationary source and sink quantities are denoted a (aj and o~b are the 
rate of formation and breakdown of the dimer in unit time and unit volume 
respectively throughout the cell. The conservation of mass equation for atoms 
and dimers read as follows, where the subscripts refer to the species label for 
the flow vector J and the concentration c: 

dcA 2 /dt = — V. Ja 2 + Of — (76 

dc A /dt = -V.J A -2a f + 2a b (30) 

The steady state condition is V. J a — ~2(ctj — 0]f) = — 2oy and V. Ja 2 — °V 
; (07 — (Jb) = o~ r where ay is a scalar flux and at thermodynamical equilibrium, 
a r = strictly. If the PLE were strictly valid then the J a, Ja 2 fluxes must 
vanish; clearly here, this is not the case. To check for flux conservation, the 
divergence term is discretized by integration over one layer, using the trapezoidal 
rule, where for any layer i, 

f V.J A2 dV = M?) = g f = 1))AV = J A2tdif (i) = J A2 (i)-J A2 (i-l) (31) 

where the layer has volume AV. Similarly, for the atomic fluxes, 

JA,dif {i) = JA{i) - J A (i - 1) = -M») + oy(i - 1))AV (32) 

leads to 

J d (i) = 2J A2 ,dif(i) + JA,dif(i) = (33) 

The plot of 3d given in 133|) in Fig. (J2U complies with the conservation law rather 
well, within statistical error. We have therefore shown that PLE is not a rigor- 
ous principle from numerical simulation where a counter-example is given, and 
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that local stochastic equilibrium dynamical variables do not necessarily have 
Gaussian (Canonical) distributions as demanded by some specialists |28l 1291 1377] 
in their theories. 



0.1 
0.08 
0.06 
0.04 
0.02 


-0.02 
-0.04 



-0.06 



J flux(Case2) 



J A flux(Case 2) 



I M I* 1*1 1 1 L I L I L ! U U U U 



"10 20 30 40 50 

X 



60 



Figure 28: Extreme thermodynamical conditions leading to the presence of 
steady state atomic and dimer fluxes. The data points are used to construct 
the difference equation in the text to verify the conservation law. 



5 Conclusion 

This study shows that the model of the molecule utilizing switching poten- 
tials does lead to typical behavior predicted from standard theories for unusual 
hysteresis-type reaction mechanisms which theorists have largely ignored, due 
perhaps to the influence of '"time-reversible" 'symmetry concepts. It is demon- 
strated that microscopic loop-like pathways does not influence the macroscopic 
thermodynamical results in any fundamental way. The method used here to 
reduce expensive 3-body calculations to easier 2-body calculations may be used 
as a basis for non-equilibrium simulation applications, which will be the subject 
of further investigations. The two body potentials yield extremely good ther- 
modynamic results whilst being super-efficient in reducing computational costs 
because the use of switches and algorithms that can preserve momentum and 
energy during potential transitions, and it is expected that semi-quantitative 
results at least can be determined for any molecular potential that is known. 
The NEWAL algorithm is effective for the extreme conditions of the simulation, 
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Figure 29: Test of divergence theorem for mass conservation via a difference 
equation. 



and would prove to be a valuable tool in reducing errors attributable to switch- 
ing potentials. A whole generation of scientific literature has been devoted to 
establishing necessary connections between the direction of material flow (mi- 
croscopic reversibility or "'time reversibility"') and thermodynamics, but the 
results here suggests that there need not be any necessary connection between 
the two. Lastly, it is shown through counter-examples that the PLE and the 
canonical averaging assumption used in recent thermodynamical theories are not 
strictly correct since internal variables do not have the flgebraic structure 

as the variables that are explicitly featured in the system Hamiltonian. We 
have demonstrated that there is a feasibility of developing an extended theory 
of equipartition (where the temperature parameters associated with any species 
motion need not be fixed and of the same value as the system thermodynamical 
temperature) on the basis of the shape of the energy distributions. It would be 
of interest to repeat and compare some of the above calculations for a conven- 
tional system without hysteresis to rule out any necessary connection between 
dynamics and equilibrium thermodynamic properties. 
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